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Solution of 2D Boussinesq systems with FreeFem++: 
the flat bottom case 

Georges Sadaka* 

o 

On 

Abstract - We consider here different family of Boussinesq systems in two 
space dimensions. These systems approximate the three-dimensional Euler equations 
and consist of three coupled nonlinear dispersive wave equations that describe 
propagation of long surface waves of small amplitude in ideal fluids over a horizontal 
bottom and which was studied in [7,9,10]. We present here a FreeFenr^ code aimed 
at solving numerically these systems where a discretization using Pi finite element 
for these systems was taken in space and a second order Runge-Kutta scheme in 
time. We give the detail of our code where we use a mesh adaptation technique. 
An optimization of the used algorithm is done and a comparison of the solution for 
different Boussinesq family is done too. The results we obtained agree with those of 
^ the literature. 

Keywords: Boussinesq systems, KdV-KdV, BBM-BBM, Bona-Smith, adaptmesh, 
lO finite element method, FreeFem++. 

o 

^— i 1. Introduction 

> 

It has often been observed that variations of the bottom could influence the damping 
of the waves including extreme ones as Tsunamis: the coral reef or the underwater 
forests in the first shoreline, mangroves; these underwater reefs are also used to 
prevent corrosion effects of coastal (see R Azerad et al. [1] and [2]). In these cases, 
the underwater relief damped the wave energy, in contrast, in other situations we 
seek to harness this energy: some companies even offer projects underwater reefs for 
erectile produce energy from waves (see 
http : / /www . aquamarinepower . com/). 
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In a continuum approximation, waves on the surface of an ideal fluid under the force 
of gravity are governed by the Euler equations. These are expected to provide a good 
model of irrotational waves on the surface of water, say, in situations where 
dissipative and surface tension effects may be safely ignored, see [16]. 

In this paper, the attention is given to a multi-dimensional Boussinesq systems which 
describes approximately the propagation of small amplitude and long wave-length 
surface waves in a three-dimensional wave tank filled with an irrotational, 
incompressible and inviscid liquid under the influence of gravity, moving and/or 
variable bottom and surface pressure. The full Boussinesq systems in 2D has been 
derived in [18]. 

Chen, Goubet, Dougalis, Mitsotakis and Saut have considered 2D models of 
Boussinesq systems with flat bottom [7,9] then with variable bottom [6,12,17], on 
the other hand Dutykh, Katsaounis and Mitsotakis have developed a code in finite 
volumes for the Boussinesq systems with variable bottom in ID ([13]) and 
Mitsotakis et al in Galerkin finite elements (using B-splines [9]). 

In order to solve numerically the Boussinesq system, we will use FreeFenr^ which 
is an open source platform to solve partial differential equations numerically, based 
on finite element methods. The FreeFem^ platform has been developed to facilitate 
teaching and basic research through prototyping. For the moment this platform is 
restricted to the numerical simulations of problems which admit a variational 
formulation. 

Thus, we develop a FreeFerr^+ code for the simulation of Boussinesq equations 
with flat bottom, we first check that the simulations provided by our numerical code 
are consistent with the results of the recent literature, including the work of 
Dougalis, Mitsotakis et Saut [9,10,11]. This establishes the adequacy of the chosen 
finite element discretization. 

The article is organized as follows: first we discretize the problem in space by using 
finite element method and in time by using an explicit second order Runge-Kutta 
scheme, then we develop all the steps of the FreeFem^ code to solve the problem 
by using the technique of mesh adaptation and at the end we present some numerical 
results where an optimization of the used algorithm is done and a comparison of the 
solution for different Boussinesq family is done too. 

2. The Problem 
2.1. Problem settings 

The 2D Boussineq system describe the surface wave propagation of small amplitude 
and large wave length. When considering an incompressible fluid flows in £1 C R 2 , 
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they are expressed as (see [16]): 

Tfc + V-V + V-(TjV)+aAV-V-fcAifc = 0; 

(2.1) 

V, + V77 + iV|V| 2 + cAVT7-dAV, = 0, 

The variables in (2.1) are non-dimensional and unsealed: X = (x,y) E £2 and t > 
are proportional to position along the channel and time, respectively, T] = 7j(X,f) is 
proportional to the deviation of the free surface from its rest position, V = V(X 5 f) = 
u(X,t) 



v(X,0 



(w, v) y = (w;v) is proportional to the horizontal velocity of the fluid at 



fd*\ /*\ 
some height, V* = ( I is the gradient, V • ( ^ 1 = d x *+d y * is the divergence 

and A* = d xx * + d yy * is the laplacian. The coefficients a, c and d are given by the 
following formulas: 

(2.2) 

where v, are real constants and ^ # ^ 1. 

We note that the dispersive constants a, b, c and d satisfy the physical constraints (see 
[3] for detail): 

a + b + c + d= ^ andc + J^O (2.3) 

We now list some of the different family of Boussinesq systems 2D in Table 1 of the 
form (2.1): 



System 


V 2 


V 


M 


References 


BBM-BBM 


2/3 








[6,7,9,10,11]. 


Bona-Smith 


2/3 < tf 2 < 1 





4-6tf 2 


[7,9,10,11]. 


3(1 -# 2 ) 


" General " Boussinesq 


< 2 < 1 


any 


any [7,9,11]. 


KdV-KdV 2/3 1 


1 


[14] 



Table 1: Examples of Boussinesq systems in 2D. 



In [9], V. Dougalis, D. Mitsotakis and J-C. Saut have studied the Well-Posedness of 
the Boussinesq systems (2.1) (where b ^ and d 7^ 0) and have shown that this 
system is at least nonlinearly well-posed locally; and in the case of KdV-KdV 
system (where b = d = 0, a = c = 1/6), F. Linares, D. Pilod and J-C. Saut proved 
recently in [14] the Well-Posedness of this system. 
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Following [19], we can write (2.1) as: 



T-AV 







Tfc + V-V + V-(TjV) + aV-r-fcATfc 







(2.4) 



0- At] 







V f + VT] + 



±V|V| 2 + cV@-dAV, 



where T= (T l J 2 ) T = (T^T 2 ). 



3. Numerical Scheme 

In this section, we present the spatial discretization using finite element method with 
Pi continuous piecewise linear functions as shown in [19] and for the time marching 
scheme an explicit second order Runge-Kutta [8] scheme as used in [9]. 
We will use in our code a mesh adaptation technic that we can use solving the 
problem by using the method based on the declaration of the problem obtained by 
the weak formulation of the system (2.4); or by using the second method that consist 
to build matrices and vectors to solve the direct system AX = B, where the matrix A 
and the vectors X,B will be defined in the sequel. 

3.1. Spatial discretization 

We let £1 be a convex, plane domain, let denote a regular, quasi uniform 
triangulation of £1 with triangles of maximum size h < 1 [5], let Vh = {vh E 
C°(£i);vh\T E Pi(r),Vr E T^} denote a finite-dimensional subspace of H 1 ^) = 
{u E L 2 (Q) s.t. ^, ^ E L 2 (£l)} where Pi is the set of polynomials of R of degrees 
< 1 and let (•; •) denote the L 2 inner product on £1. 

Consider the weak formulation of the system (2.4), find r\h^ u h^h E Vh such that 
V% E Vh we have: 

(Tj; (p h ) - (Au h \ (p h ) = 0; (T 2 ; (p h ) - (Av h ; q> h ) = 0; (& h ; (p h ) - (Arj^; (p h ) = 0: 
({Id - bA)riht + V • (u h ;v h ) + ri hx u h + r\ h u hx + r\ hy u h + r\ h u hy + aV • (T^;T|) ; <p h \ 0: 

( (Id - dA)u ht + Tfhx + u h Uhx + v h v hx + c©^; (p h ) = 0: 



(Id-dA)v ht + ri h y + u h u h y + v h v h y + c® h y\(p h ) = 0. 



(3.1) 
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To simplify, we denote <I> = <p h , E = r\ h , U = u h , V = v h , T = & h , P = and Q = Tj r , 
so the system (3.1) is equivalent to the following system: 

(P;<D) = (AU;<D); (Q;<D) = (AV;<D); (T;<D) = (AE;<D); 
(Id-bA)d t E;^ = -(V-(U;V)+E x U + EU x + E y V + EV y + aV-(P;Q);$) 

- -F(E,U,V,P,Q); 
(Id-dA)d t \J;<£>\ = -(E x + UU x + VV x + cT x ;<£) = -G(E,U,V,T); 

(Id-dA)d t Y;<t>) = -(E y + UU v + VV y + cT y ;<£) = -H(E,U,V,T). 

(3.2) 

3.2. Time marching scheme 

Our method is based on an explicit second order Runge-Kutta scheme. To this end, let 
us denote by (E" +1 ,U"+\ V" +1 ) and (E>\U", V",P n ,Q",T n ) the approximate value 
at time t = t n+l and / = t n , respectively and by St the time step size. Then, by using 
(3.2), the unknown fields at time t = t n+l are defined as the solution of the system 

(F»;4>) = (AU";<I>) ; (Q";<I>) = (AV";<I>) ; (T";d>) = (AE";<&) ; 
(E" +1 ;0) = (E" + 1 ;0);(U" +1 ;0) = ^ U ;<£); 

(yn+l.fy = (yn + Z .^y 

(3.3) 

where: 

' ((W-M)E w ;0) = -Sr-F(E",U",V n ,P n ,Q"); 
' (Id - dA)U kl ; <£) = - Sf • G (E" , U" , V" , T" ) ; 
(Id - dA) V a ; <J>) = - Sr • H (E» , W , V" , T" ) ; 
[ (P* 1 ;^) = <U» +U«;0>;(Q«;<I>> = <V« +V yy ;<I>> ;(T«;0> = <E» +E«;d>>. 

(3.4) 

and 

((W-&A)E°;0) = -5f-F(E" + E w ,U" + U H ,V" + V H ,P" + P w ,Q" + Q w ); 
((W-dA)U* 2 ;<I>) = -5?-G(E"+E M ,U n + U fcl ,V" + V fcl ,T n + T fcl ) 
((/</- </A)V H ;<I>) = -St-n(E n + E kl ,V n +U H ,\ n + V H ,T n + T kv 

(3.5) 

By integrating by parts where we have second order derivative and by developing all 
the terms of the first order derivative in (3.3), (3.4) and (3.5), we deduce: 

/ d\J" \ I dV n 

(P»;0) = - (VU«; V<D) + < — ;<!> ) ; (Q»;4>> = - (W; V4>) + ( ^— ;<I> 

\ on I da \ on i da 



' dE n 

(T n ;4>) = - (VE B ; V<J>) + ( — ;<!> 

' an 

(3.6) 



<9n 
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da 



da 



(E kl ;<t>)+b(VE kl ;V<t>)-b{-j—;<i> 

(U kl ;<Z>)+d ( VU* 1 ; VO> - d / ^ 

(V* 1 ; O) + J ( VV* 1 ; VO) - </ / 

/<9U a 

(p«;<I>) = _(vU a ;V<I>) + / -j- ;0 

/<9V a 

(Q";0) = - (VV^ 1 ; VO) + ( — — ;<!> 

\ c,n / 
' <9E 



an 



(T H ;0> = - (VE* 1 ; V<J>> + ( 



da 



-5f-F(E",U",V",P",Q' 1 ); 

-5f-G(E",U",V",T"); 

-5f-H(E",U",V",T"); 



and 



(3.7) 



/ (9E^ 2 

(E* 2 ; + i < VE* 2 ; VO> - Z> ( -3— ; O ' 

\ dn I da 

-5f-F(E" + E H ,U" + U H ,V" + V H ,P" + P H ,Q" + Q H ); 
/ d\J k2 \ 

(\J k2 ;<i>)+d{VU k2 ;V<i>)-d(^—;<i>) = 

\ ° n /da 
-5f-G(E" + E H ,U' , + U a ,V" + V a ,T" + T /:1 ); 

/ dV k2 

(V k2 ;<t>)+d(Vy k2 ;V<i>)-d( ;0 

\ dn / da 
-5f-H(E" + E H ,U" + U a ,V n + V a ,T' ! + T a ). 

(3.8) 

Remark: It's easy with FreeFenr^ to define boundary condition, in fact if we have 
the Dirichlet Boundary Conditions on a border Ti C M like = /, then it is defined 
as on ( gamma 1 , u=f ) , where u is the unknown function in the problem. We note 

dU, 

that the Neumann Boundary Conditions on T2 C R, like r? = g> appear in the 

dn 

Weak formulation of the problem after integrating by parts for example in the system 

(3.6) we have ( — ; <I> \ = (g; <&) T2 = / g • <& which is defined in FreeFeiTH-+ by 

\on / y 2 2 Jr 2 

intld (Th, gamma2 ) (g*phi) where Th is the triangulated domain of £1. We 
will see in the next section how it's also easy to define the Bi-Periodic Boundary 
Conditions. 
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We remark also that the system (3.3) can be written on the following matrix form: 



M \ /E" +1 \ 
M • U" +1 
M J \ V" +1 / 



(E" + 
(U" + 
(V" + 





+ E k2 




2 




2 




+ y *2 



;<*>> 



(3.9) 



where = / <pi<p/djcdy is the mass matrix. 



Algorithm 1: Finally, to solve the systems (3.3), (3.6), (3.7) and (3.8), we follow as: 

E n = E° = rfoo = rio,V n = U° = W/z0 = no, V" = V° = v h0 = v 
P" = flw, Q* = e*o, T 1 = 7m), P* 1 = P M1 , Q* 1 = G«i , T* 1 = T hkl 
E n+1 = ri h ,V n+1 =u h ,\ n+1 =v h 



Set 
Set 
Set 
Set 

For t = 0:St:T 



Set 



Mesh adaptation, (optional) 

Compute Pm.QmJm Compute rf hku u hkh v hk i 

Compute P^kuQhkuThki Compute T] M2 ,wm2,Vm2 

Compute Tl hl u h ,v h 

rim = rih,u h0 = u h ,v h0 = v h 



End for 



Algorithm 2: Another method to solve the systems (3.6), (3.7) and (3.8), taking into 
account (3.9): 



Set E n = E° = rfoo = T]o,U" = U° = u h0 = u ,Y n = V° = v h0 = v 

Set P* = P h0 ,Q n = G«>, = r, , P M = P hk i , Q M = Q hk i , T* 1 = r M1 

Set E n+l = rj^,U n+1 = u h ,\ n+l = v h 

Set E* 1 = ffo^U* 1 = M/wki, V w = v hkh E k2 = r] M2 ,U^ = w M 2, V* 2 = v hk2 

Compute A (if we want to use the mesh adaptation 

we must compute A in the for-loop time) 



For t = 0:8t:T 

Mesh adaptation, (optional) 

Update rjfco t?m); w M) ^m'^m vw'^h Vh'^h W^h — v h\ 
(with mesh adaptation) 

Compute Pm,Qm,Tm Compute ri hku u hku v hk i 
Compute PticuQhkuThki Compute rf hk2 ,u hk2 ,v hk 2 
Compute A (with mesh adaptation) 
Set X = [r\hi u h, Vh] Compute B Solve AX = B 
Set rfho ri h ,u h0 u h ,Vho Vh 
End for 
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4. Code 

In this section we will present by details all the step of the FreeFem^ code to solve 
(3.3) to (3.9). 

4.1. Declaration of the problems 

Note that in FreeFem^ the scalar product inL 2 : (., %) = / . • q> h = int 2 d ( Th ) ( 

. *phih ) ; also we can define a macro for the right hand side function F (E, U, V, P, Q), 
G(E,U, V,T) ,H(E,U, V,T) defined in (3.2) using the keyword macro, that will 
be used in the sequence, as: 

macro F(e,u,v,p,q) (div (u, v) +dx (e ) *u+e*dx (u) +dy (e ) *v+e*dy ( v) + 

**a*div (p, q) ) // 
macro G(e,u,v,t) (dx (e ) +dx (u) * (u) +dx ( v) * ( v) +c*dx ( t ) ) / / 
macro H(e,u,v,t) (dy (e ) +dy (u) * (u) +dy ( v) * ( v) +c*dy ( t ) ) / / 



We note that all the variable (e , u , v, p , q, t) used in the macro are dummies. 
We declare the problem for XJ kl defined is the system (3.7) and for E" +1 defined in 
the system (3.3) as: 

problem UHK1 (uhkl , phih) = int2d(Th) (uhkl*phih) + int2d(Th) ( 
**grad (uhkl) ' *grad (phih) *d) + int2d(Th)( G (etahO , uhO , vhO , 
^ThO) *phih*dt ) +"Boundary Conditions of uh for uhkl"; 

problem ETAH (etah, phih) = int2d(Th) (etah*phih) - int2d(Th) ( 
*>etahO*phih) - int2d (Th) ( (etahkl + etahk2) * phih 12.)} 



The declaration of the problems for all other variables are written in the same form. 

Remark: In order to make our code faster, we can use the keyword in it in the 
declaration of the problem. When init=0 the mass matrix is computed and when 
init=l the mass matrix is reused so it is much faster after the first iteration. 

4.2. Solve of the problems 

To solve all the problems defined above, we make a for-loop time and we call the 
problems by their names when we want them to be solved, then we update the data 
and at the end we plot the solution using the keyword plot. 
We note that in each iteration of the for-loop a mesh adaptation will be done which 
depend on the error (err) which is the Pi interpolation error level, where hmin is 
the minimum edge size and nbvx is the maximum number of vertices generated by 
the mesh generator. 



for (real t=0 . ; t<=T; t+=dt) { 
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} 



Th=adaptmesh (Th, etahO , err=le-4 , hmin=Dx, nbvx=le6 ) ; // we 

^can use adaptmesh each 10 iterations or more. 
PHO; QHO; THO; 

ETAHK1; UHK1; VHK1; 

PHK1; QHK1; THK1; 

etahOpkl=etahO+etahkl; uhOpkl=uhO+uhkl ; vh0pkl=vh0+ 
**vhkl; 

Th0pkl=Th0+Thkl; PhOpkl=PhO+Phkl ; Qh0pkl=Qh0+ 

**Qhkl; 

ETAHK2; UHK2 ; VHK2 ; 

ETAH; UH; VH; 

etahO=etah; uh0=uh; vh0=vh; //update of the data 

plot (etahO , cmm=" t = "+t + " sec" , f ill=true, value=true, dim=3 ) ; 



In order to use the second method, we build the matrix A before the for-loop time as: 

varf Mass(u, phih) = int2d(Th) ( u * phih ); 
matrix A, MASS; 
MASS = Mass (Vh, Vh) ; 

A= [[MASS, 0, 0],[0, MASS, 0],[0, 0, MASS]]; 
set (A, solver=GMRES ) ; // to be set 

Then we build the vector B in the for-loop time as: 

for (real t=0 . ; t<=T; t+=dt) { 
PHO; QHO; THO; 

ETAHK1; UHK1; VHK1; 

PHK1; QHK1; THK1; 

etahOpkl=etahO+etahkl; uh0pkl=uh0+uhkl ; vhOpkl=vhO+vhkl ; 
Th0pkl=Th0+Thkl; Ph0pkl=Ph0+Phkl ; QhOpkl=QhO+Qhkl ; 

ETAHK2; UHK2 ; VHK2 ; 

Vh Bl, B2, B3, etahklpk2D2, uhklpk2D2, vhklpk2D2; 
real[int] B ( 3^Vh . ndof ) , X ( 3*Vh . ndof ) , X0 ( 3 *Vh . ndof ) , W(3^Vh. 
^►ndof ) ; 

etahklpk2D2 = . 5*etahkl + .5*etahk2; 

uhklpk2D2 = .5^uhkl + .5*uhk2; 

vhklpk2D2 = .5^vhkl + .5*vhk2; 

X0=[etah0[], uhO[] A vhO [ ] ] ; 

Bl [ ] =MASS*etahklpk2D2 [ ] ; 

B2 [ ] =MASS*uhklpk2D2 [ ] ; 

B3 [ ] =MASS*vhklpk2D2 [ ] ; 

B=[Bl[] f B2[] f B3[]]; 

X = A^-l^B; 

W = X + X0; 

[etah[] , uh[], vh [ ] ] = W; 

etahO=etah; uh0=uh; vh0=vh; //update of the data 

plot (etahO, cmm=" t=" +t+" sec" r f ill=true A value=true, dim=3) ; 
} 

Finally, if we want to use mesh adaptation in the second method, we must compute 
the matrix A in the for-loop time. 
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4.3. Numerical simulations 

In the sequel, we present the results of numerical simulations of the evolution of 
initially localized heaps of fluid of initial velocity zero. Unless specified, all com- 
putations were performed on the square £1 =] — 40,40[x] — 40,40 [, a Pi continuous 
piecewise linear functions was used for the finite element space and for all the nu- 
merical simulations, we work with the space discretization Ax = 0.5 and the time 
step At = 0.1. 



4.3.1. Rate of convergence. At the beginning, we prove in the figure below, that 
the RK2 time scheme considered for the BBM-BBM system is of order 2. In this 
example, we took zero Dirichlet homogenous Boundary Conditions for T]^, Uh and 
on the whole boundary and we have consider the following exacts solutions: 

r] ex = ^ • sin(nx) - (y - I) -y, 

u ex = e* - x-cos(37rx/2) • sin(7ry), 

v ex = e* -sm(7tx) • cos(37ry/2) -y. 

Then, we compute the corresponding right hand side in order to obtain the L 2 norm 
of the error between the exact solution and the numerical one in the table below. 



N 


\f]h — ^lex | L 2 


\Uh — "ex L 2 


\Vh Vex | L 2 


10 


0.00871494 


0.0233966 


0.0230945 


20 


0.00265707 


0.00641675 


0.00632314 


40 


0.000670301 


0.00160223 


0.00157848 


80 


0.0001817 


0.000419198 


0.000412791 


160 


4.80657e-05 


0.000108456 


0.000106767 



Table 2: L 2 norm of the error for T] , w, v. 



4.3.2. Computation time. In this section, we consider the BBM-BBM Boussinesq 
system on the domain £1 =] — 40, 40 [ x ] — 40, 40 [ with Pi continuous piecewise linear 
functions, the space discretization Ax = 0.5, the time step At = 0.1 and as initial data 
^hoi^^y) = 0.2e~( x2+y2 ^ 5 , Uho(x,y) = Vho(x,y) = with zero Dirichlet homogenous 
Boundary Conditions for 7]/^, Uh and Vh on the whole boundary. 
In order to solve this system, we will show the time comparison of different method: 

• Ml to solve Algorithm 1 without using adaptmesh technique and the keyword 

init. 
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i loglog scale tor BBM-BBM Bot 



10 ^L , , , , , , , 

10 10 

Step time dt 

Figure 1: Rate of convergence for BBM-BBM. 



• Mlinit to solve Algorithm 1 using the keyword in it and without using 
adaptmesh technique. 

• M1A-4 to solve Algorithm 1 using adaptmesh technique with err=le-4 
and without the keyword init. 

• M1A-2 to solve Algorithm 1 using adaptmesh technique with err=le-2 
and without the keyword init. 

• M2 to solve Algorithm 2 without using adaptmesh technique and the keyword 

init. 

• M2init to solve Algorithm 2 using the keyword init and without using 
adaptmesh technique. 

• M2A-4 to solve Algorithm 2 using adaptmesh technique with err=le-4 
and without the keyword init. 

• M2A-2 to solve Algorithm 2 using adaptmesh technique with err=le-2 
and without the keyword init. 



We present in Table 3, the time of computation in second at time T = 10s using all 
the different method cited before. All computation was made on a Macbook OS X, 
Intel core 2 Duo (CPU), 4Go (Memory), 2 Ghz (Processor). 



We note that, without using the mesh adaptation technique, we have the same result 
for all the computed solution, so we can see from Table 3 that the best method to use 
is the M2init. 
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Ml 


Mlinit 


M1A-4 


M1A-2 


1555.48 


722.507 


115.485 


45.3462 


M2 


M2init 


M2A-4 


M2A-2 


1217.01 


619.462 


99.5689 


39.3365 



Table 3: Comparison of computation time for the different method used to solve the 
BBM-BBM system. 

In other hand, using the mesh adaptation technique, we can remark from Table 3 that 
the computation time is better then other method, but unfortunately, we have a little 
difference between the computed solution, that we plot in Figure 2 the square of L 2 
norm of the error between the computed solution using the Ml method and the one 
computed with the M1A-8, M1A-6, M1A-4, M1A-2 methods where we have 
err=le-8, err=le-6, err=le-4, err=le-2, respectively vs the time till 
T = 30s. We also plot in Figure 3 the mean of all the square of L 2 norm of the error 
computed for different mesh adaptation method. 



Comparison of the error between the solution using adaptmeah and without using adapt mesh 




Figure 2: Comparison of the error between the solutions. 



We can remark from this result that we have the same result using M1A-8, M1A-6, 
Ml A-4 method and the mean error between the solution computed with these method 
and the computed one using the Ml is of order 10~ 5 and we can see the large time 
difference. 
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Mean of the error between the solution using adapt mesh and without using adapt mosh 




Figure 3: Mean of the error between the solutions. 



4.3.3. Reflection of expanding symmetric waves at two boundaries of the BBM- 
BBM Boussinesq system. In Figures 4 and 5 we show the reflection from two 
parts of the boundary of an expanding symmetric wave of the BBM-BBM Boussinesq 
system where a = c = and b = d = 1/6. For this experiment we used as initial 

data the functions T}ho(x,y) = .2e~( x2+y2 ^ 5 , Uho(x,y) = Vho(x,y) = 0. We used zero 
Neumann Boundary Conditions for 7]^ on the whole boundary, zero Dirichlet data 
for Uh and Vh on x = —40 and y = 40 (where we have the wall), and zero Neumann 
boundary data for Uh and v^oni^ 40 and y = —40. The expanding waves are 
reflected from the x = — 40 and y = 40 parts of the boundary. 
We note that in Figure 4 we show the effect of the mesh adaptation following the 
evolution of r\h in time and in Figure 5 we show the propagation of the solution T]^. 
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■■■■ 



_________ 



Figure 4: Propagation of the mesh for the BBM-BBM Boussinesq system where 
a = c = and b = d = 1 /6 for different time r = {0.1,20, 40, 70} 




Figure 5: Propagation of the solution of the BBM-BBM Boussinesq system where 
a = c = and b = d = 1 /6 for different time t = {0, 20, 40, 70} 
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4.3.4. Expanding symmetric waves under the KdV-KdV Boussinesq system. 

In Figure 6, we present the evolution of the r\h profile emanating from the radially 
symmetric initial data riho(x,y) = .5e~( x2+y2 ^ 5 , itho(x,y) = Vho(x,y) = 0, under the 
KdV-KdV Boussinesq system where a = c = 1 /6 and b = d = 0. We used Bi-Periodic 
Boundary Conditions for T]^, Uh and Vh and we work with the time step At = 0.001. 
We remark here that with these Bi-Periodic Boundary Conditions for 7], u and v and 
their derivatives, in addition by integrating the equations in the system (2.1) on the 

hole domaine, we deduce the following mass conservation: (Id — bA) / rf t = and 

Jq 

the relations (Id — dA) / u t = 0, (Id — dA) I v t = 0. Hence: 
Jo. Jq. 

/ 7]=cte= / Tjo, / ii = etc = / U(). / v — cte— I vq. (4.1) 
Jq Jq Jq Jq Jq Jq 

In other hand, numerically, we see that these defined quantity are well conserved 
over time and we have: 

/ 7] = etc = / 77o = 7.84527, / u = cte= w = 0, / v = cte= v () = 0. 
Jq Jq Jq Jq Jq Jq 





Figure 6: Propagation of the solution of the KdV-KdV Boussinesq system where 
a = c = 1/6 and b = d = for different time t = {0, 10, 20, 60} 



We can see in Figure 6 from t = 10 a small amplitude periodic profile (ripples) 
which are propagating in front of the wavefront and which has been observed in [4] 
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(in the case of ID KdV-KdV system). These ripples are still observed when using P2 
elements but they do not infer to the numerical stability, we still have the mass 
conservation of the fluid. 

Other simulation for different Boussinesq systems can be found in [18]. 



4.4. Comparison of the results with Mitsotakis et al. 

We validate also our code by comparing the results with those of Mitsotakis et al. 
[9] (in Figure 8, page 847) where they consider the cross section in x direction to 
the x\ component of the solution from radially symmetric initial data of the form 
7] h0 (x,y) = .2e~^ x2+y2 ^ 5 with u h0 (x,y) = v m {x,y) = under the BBM-BBM and the 
Bona-Smith system, both considered with zero Dirichlet boundary conditions for 
r/ , u and v. We remark in Figure 7, that the shape of the solution for the results of our 
code at the left part and for those of Mitsotakis et al. code at the right part are similar. 
We note also that in their paper, they used Galerkin finite elements based on tensor 
products of smooth splines and an explicit second order (for BBM-BBM system 
with bilinear splines) and fourth order (for Bona-Smith system with bi-cubic splines) 
Runge-Kutta scheme while in the present work, we used Pi continuous piecewise 
linear functions and an explicit second order Runge-Kutta scheme for all systems. 



4.5. Comparison of different Boussinesq models 

We compare here KdV-KdV, BBM-BBM and Bona-Smith models as defined in 
section 2.1 Table 1. 

In Figures 8, we present a comparison of the evolution of the T] component of the 
solution from radially symmetric initial data of the form r\m( x ->y) = -5e~^ x +y ^ 5 
with Uho(x,y) = v/jo(x,y) = 0, under the BBM-BBM (solid line with zero Dirichlet 
b.c, At = 0.1, Ax = 0.5), the Bona-Smith system (dotted line with zero Dirichlet 
b.c, At = 0.1, Ax = 0.5) and the KdV-KdV system (dashed line with periodic b.c, 
At = 0.001, Ax = 0.5). Figure 8 shows the cross sections of the T] profiles for 
different time in the x - direction. The speed and the amplitude of the outgoing 
front is approximately the same for the BBM-BBM and Bona-Smith systems but 
the pattern of the oscillations behind the fronts are different: in the case of the Bona- 
Smith system the two outgoing wave trains have practically separated by t = 25s, 
while the larger in amplitude dispersive oscillatory tails of the BBM-BBM solution 
seem to be still interacting. In other hand, the speed of the outgoing front for the 
KdV-KdV system is approximately the same with other systems while the amplitude 
for the internal crest is bigger, we remark also that the solution still interacting after 
t 25s. 
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Figure 7: Cross sections in x direction of r\{x,y,t) at / = 30s, Bona-Smith (above) 
and BBM-BBM (below), using our code (left) and Mitsotakis et al. code (right) 
borrowed from the article [9]. 



5. Conclusion 

We have presented a numerical approach with FreeFenr^ to solve the Boussinesq 
systems with a flat bottom, we validated our code and establishes the adequacy of the 
chosen finite element discretization by comparing the results with those of Mitsotakis 
et al. We have established also the feasibility of simulating complex equations of 
hydrodynamics as Boussinesq systems with FreeFeiT^ and we have optimized the 
algorithm that we use in section 4.3.2. 

Using this approach, we can consider the case of a variable bottom (in space and/or in 
time), see [18] which is an ongoing work to appear soon. As a feature we address the 
simulation of Tsunamis with our approach, by including realistic data (bathymetry, 
generation of tsunami waves). 
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Figure 8: Comparison of the cross sections in x direction of t]{x,y,t) at time / = 
{2, 5, 10, 15, 20, 25}s, where in all the figure, the dotted line is for Bona-Smith system, 
the dashed for KdV-KdV one and the solid line for BBM-BBM one. 



Acknowledgements: This work is supported by the regional program "Appui a 
lemergence" of the Region Picardie. I would like to thank my PhD advisor Jean- 
Paul Chehab (LAMFA, Amiens), Denys Dutykh (LAMA, Savoie), Gloria Faccanoni 
(IMATH, Toulon), Kirill Pichon Gostaf (LJLL, Paris), Frederic Hecht (LJLL, Paris), 
Antoine Le Hyaric (LJLL, Paris), Youcef Mammeri (LAMFA, Amiens), Olivier 
Pantz (CMAP, Paris) and Dimitrios Mitsotakis (IMA, University of Minnesota) for 
fruitful discussions and remarks. 



Solution of 2D Boussinesq systems with FreeFem++: the flat bottom case19 



References 

1. Pascal Azerad, Frederic Bouchette, Benjamin Ivorra, Damien Isebe and 
MOHAMMADI. Shape optimization of geotextile tubes for sandy beach protection. Int. J. Numer. 
Meth. Engng, Vol. 74, pp. 1262-1277, 2008. 

2. Pascal Azerad, Frederic Bouchette, Damien Isebe and Bijan Mohammadi. 
Optimal shape design of defense structures for minimizing short wave impact. Coastal 
Engineering, Vol. 55, pp. 35-46, 2008. 

3. Jerry Lioyd Bona, Min Chen and Jean-Claude Saut. Boussinesq equations and other 
systems for small amplitude long waves in nonlinear dispersive media: II. the nonlinear theory. 
Nonlinearity, 17, no3, 925-952, 2004. 

4. J. L. Bona, V. A. Dougalis and D. E. Mitsotakis. Numerical solution of KdV-KdV 
systems of Boussinesq equations I: The numerical scheme and generalized solitary waves. Math. 
Comput. Simulation, 74:214-228, 2007. 

5. Suzanne C. Brenner and L. Ridgway Scott. The Mathematical Theory of Finite Element 
Methods. Springer-Verlag, New York, 1994. 

6. MlN Chen. Numerical investigation of a two-dimensional Boussinesq system. Discrete and 
Continuous Dynamical Systems, Vol 23, no4, 1169-1190, April 2009. 

7. Min Chen and Olivier Goubet. Long-time asymptotic behavior of two-dimensional 
dissipative Boussinesq systems. Discrete and Continuous Dynamical Systems Series S, Vol 2, 
nol, 37-53, March 2009. 

8. Jean-Pierre Demailly. Analyse numerique et equations differentielles. Presses Uni- 
versitaires de grenoble, 1991. 

9. Vassilios Dougalis, Dimitrios Mitsotakis and Jean-Claude Saut. On some 
Boussinesq systems in two space dimensions: theory and numerical analysis. M2AN, no. 5, 
825-854, Vol 41, 2007. 

10. Vassilios Dougalis, Dimitrios Mitsotakis and Jean-Claude Saut. On initial- 
boundary value problems for a Boussinesq system of BBM-BBM type in a plane domain. AIMS' 
Journal, Vol 23, 2009. 

11. Vassilios Dougalis, Dimitrios Mitsotakis and Jean-Claude Saut. Boussinesq 
systems of Bona-Smith type on plane domain: Theory and Numerical Analysis. Journal of 
Scientific Computing, Vol. 44, no2, pp. 109-135, 2010. 

12. Denys Dutykh and Frederic Dias. Dissipative Boussinesq equations. C. R. Mecanique, 
335,559-583,2007. 

13. Denys Dutykh, Theodoros Katsaounis and Dimitrios Mitsotakis. Finite volume 
schemes for dispersive wave propagation and runup. Computational Physics, 230 (8), 3035 - 
3061,2011. 



20 



Georges Sadaka 



14. Felipe Linares, Didier Pilod and Jean- Claude Saut. Well-posedness of strongly 
dispersive two-dimensional surface waves Boussinesq. ArXiv: 1 103. 41 59v2, 1 1 Apr 201 1. 

15. Brigitte Lucquin and Olivier Pironneau. Introduction to Scientific Computing. Wiley, 
1998. 

16. Alain Miranville and Roger Temam. Mathematical modeling in continuum mechanics. 
Cambridge University Press, 2005. 

17. DIMITRIOS Mitsotakis. Boussinesq systems in two space dimensions over a variable bottom 
for the generation and propagation of tsunami waves. Mat. Comp. Simul., 80:860-873, 2009. 

18. Georges Sadaka. Etude mathematique et numerique d' equations d'ondes aquatiques amorties. 
These de VUniversite de Picardie Jules Verne - Amiens, 201 1. 

19. Mark Walkley and Martin Berzins. A finite element method for the two-dimensional 
extended Boussinesq equations. Int. J. Numer. Meth. Fluids, 39:865-885, 2002. 



